

summary2 <- function(model, clust = "clust3"){
  modses <- sqrt(diag(vcov(model))) 
  modcoefs <- summary(model)$coefficients; modcoefs[,2] <- modses; 
  modcoefs[,3] <- modcoefs[,1]/modcoefs[,2]; modcoefs[,4] <- 1 - pnorm(abs(modcoefs[,3]))
  return(modcoefs)
}

myround <- function(vec, dig){return(formatC(round(as.numeric(as.character(vec)),dig),dig,format="f"))}

# for writing tables
ptable <- function(contents, file){
  print(xtable(contents), type='latex', sanitize.text.function=identity, include.rownames = FALSE, include.colnames = FALSE, 
  file = file, only.contents = TRUE, hline.after = NULL, add.to.row = addtorow)
}

regstars <- function(vec){ 
  sts <- ""; if(0 < quantile(vec,  .025) | 0 > quantile(vec,  .95)){sts <- "^{*}"}
  if(0 < quantile(vec,  .005) | 0 > quantile(vec,  .975)){sts <- "^{**}"}
  if(0 < quantile(vec,  .0005) | 0 > quantile(vec,  .995)){sts <- "^{***}"}; return(sts)
}

# functions for first differences
multlogcs2 <- function(model, cur, csvar, csquant){ 
  set.seed(12345);
  if(csvar != "diff") diff <- "Differentiated"
  factvar <- class(cur[,csvar]) == "factor"
  modobj <- mod <- zelig(model, data = cur, model = "mlogit", cite = F)
  betas <- mvrnorm(10000, mu = coef(mod), Sigma = vcov(mod)); 
  sup <- c((1:(ncol(betas)/3))*3 -2); div <- sup + 1; opp <- sup + 2
  modeln <- model; lhs(modeln) <- quote(korrelus);
  mod <- lm(modeln, data = cur)
  cs0 <- ifelse(factvar, csquant[1], quantile(cur[,csvar], csquant[1])); cs1 <- ifelse(factvar, csquant[2], quantile(cur[,csvar], csquant[2]));
  ndata <- subset(cur, select = all.vars(model)[-c(1,c(1:length(all.vars(model)))[all.vars(model)%in% c("diff","country","country2")])]); names.ndata <- names(ndata)
  if(is.null(dim(ndata))){ndata <- as.numeric(ndata)}else{ndata <- apply(ndata, 2, as.numeric)}
  if(is.null(dim(ndata))){cov.values <- median(ndata)}else{cov.values <- apply(ndata, 2, median)}
  cov.values <- data.frame(t(cov.values)); names(cov.values) <- names.ndata; cov.values$diff = diff; cov.values$country <- "Korea"; cov.values$country2 <- "Korea"
  cov.values[,csvar] <- cs0; cov.values0 <- cov.values; 
  cov.values[,csvar] <- cs1; cov.values1 <- cov.values;
  covs0 <- model.matrix(mod, data = data.frame(korrelus = 1, cov.values0))
  nop.pred0 <- 0; opp.pred0 <- betas[,opp]%*%t(covs0); sup.pred0 <- betas[,sup]%*%t(covs0); div.pred0 <- betas[,div]%*%t(covs0); 
  sum.pred0 <- exp(nop.pred0) + exp(opp.pred0) + exp(div.pred0) + exp(sup.pred0);
  covs1 <- model.matrix(mod, data = data.frame(korrelus = 1, cov.values1))
  nop.pred1 <- 0; opp.pred1 <- betas[,opp]%*%t(covs1); sup.pred1 <- betas[,sup]%*%t(covs1); div.pred1 <- betas[,div]%*%t(covs1); 
  sum.pred1 <- exp(nop.pred1) + exp(opp.pred1) + exp(div.pred1) + exp(sup.pred1);
  csn <- exp(nop.pred1)/sum.pred1 - exp(nop.pred0)/sum.pred0; csnout <- paste(myround(100*quantile(csn, .5), 2), regstars(csn), sep = "")
  css <- exp(opp.pred1)/sum.pred1 - exp(opp.pred0)/sum.pred0; cssout <- paste(myround(100*quantile(css, .5), 2), regstars(css), sep = "")
  csd <- exp(div.pred1)/sum.pred1 - exp(div.pred0)/sum.pred0; csdout <- paste(myround(100*quantile(csd, .5), 2), regstars(csd), sep = "")
  cso <- exp(sup.pred1)/sum.pred1 - exp(sup.pred0)/sum.pred0; csoout <- paste(myround(100*quantile(cso, .5), 2), regstars(cso), sep = "")
  cs <- c(csoout, csdout, cssout, csnout)
  return(list(modobj, cs))
}

logcs2 <- function(model, cur, csvar, csquant, setvar = data.frame()){ 
  set.seed(12345);
  if(csvar != "diff" & ("diff" %in% names(setvar) == FALSE)) diff <- "Differentiated"
  factvar <- class(cur[,csvar]) == "factor"
  modobj <- mod <- zelig(model, data = cur, model = "logit", cite = F); 
  betas <- mvrnorm(10000, mu = coef(mod), Sigma = vcov(mod)); 
  modeln <- model; lhs(modeln) <- quote(korrelus);
  mod <- lm(modeln, data = cur)
  cs0 <- ifelse(factvar, csquant[1], quantile(cur[,csvar], csquant[1])); cs1 <- ifelse(factvar, csquant[2], quantile(cur[,csvar], csquant[2]));
  ndata <- apply(cur[,all.vars(model)[-c(1,c(1:length(all.vars(model)))[all.vars(model) %in% c(names(setvar), "diff","country","country2")])]], 2, as.numeric)
  cov.values <- apply(ndata, 2, median); cov.values <- data.frame(t(cov.values)); 
  if(ncol(setvar) > 0) {cov.values$set <- setvar[1,]; names(cov.values)[(ncol(cov.values)-ncol(setvar)+1):ncol(cov.values)] <- names(setvar)}
  if("diff" %in% names(setvar) == FALSE) {cov.values$diff = diff}; cov.values$country = "Korea"; cov.values$country2 = "Korea"
  cov.values[,csvar] <- cs0; cov.values0 <- cov.values; 
  cov.values[,csvar] <- cs1; cov.values1 <- cov.values;
  covs0 <- model.matrix(mod, data = data.frame(korrelus = 1, cov.values0)); pred0 <- betas%*%t(covs0); 
  covs1 <- model.matrix(mod, data = data.frame(korrelus = 1, cov.values1)); pred1 <- betas%*%t(covs1);
  cs <- inv.logit(pred1) - inv.logit(pred0); csout <- paste(myround(quantile(100*cs, .5), 2), regstars(cs), sep = "")
  cs <- c(csout)
  return(list(cs, quantile(inv.logit(pred0), c(.5, .025, .975)), quantile(inv.logit(pred1), c(.5, .025, .975))))
}

lmcs2 <- function(model, cur, csvar, csquant, setvar = data.frame()){ 
  set.seed(12345);
  if(csvar != "diff" & ("diff" %in% names(setvar) == FALSE)) diff <- "Differentiated"
  factvar <- class(cur[,csvar]) == "factor"
  modobj <- mod <- zelig(model, data = cur, model = "ls", cite = F); summary(mod)
  betas <- mvrnorm(10000, mu = coef(mod), Sigma = vcov(mod)); 
  modeln <- model; lhs(modeln) <- quote(korrelus);
  mod <- lm(modeln, data = cur)
  cs0 <- ifelse(factvar, csquant[1], quantile(cur[,csvar], csquant[1])); cs1 <- ifelse(factvar, csquant[2], quantile(cur[,csvar], csquant[2]));
  ndata <- apply(cur[,all.vars(model)[-c(1,c(1:length(all.vars(model)))[all.vars(model) %in% c(names(setvar), "diff","country","country2")])]], 2, as.numeric)
  cov.values <- apply(ndata, 2, median); cov.values <- data.frame(t(cov.values)); 
  if(ncol(setvar) > 0) {cov.values$set <- setvar[1,]; names(cov.values)[(ncol(cov.values)-ncol(setvar)+1):ncol(cov.values)] <- names(setvar)}
  if("diff" %in% names(setvar) == FALSE) {cov.values$diff = diff}; cov.values$country = "Korea"; cov.values$country2 = "Korea"
  cov.values[,csvar] <- cs0; cov.values0 <- cov.values; 
  cov.values[,csvar] <- cs1; cov.values1 <- cov.values;
  covs0 <- model.matrix(mod, data = data.frame(korrelus = 1, cov.values0)); pred0 <- betas%*%t(covs0); 
  covs1 <- model.matrix(mod, data = data.frame(korrelus = 1, cov.values1)); pred1 <- betas%*%t(covs1);
  cs <- pred1 - pred0; csout <- paste(myround(quantile(100*cs, .5), 2), regstars(cs), sep = "")
  cs <- c(csout)
  return(list(cs, quantile(pred0, c(.5, .025, .975)), quantile(pred1, c(.5, .025, .975))))
}
